R Notebook

sce = readr::read_rds("~/results/sub_SZ_AstOliSST.rds")
sce = as(sce, "SingleCellExperiment")

pb.logcounts

class: SingleCellExperiment 
dim: 17658 148 
metadata(3): experiment_info agg_pars n_cells
assays(7): 2 6 ... 25 30
rownames(17658): LINC00115 FAM41C ... AC004556.1 AC240274.1
rowData names(0):
colnames(148): 118847 118848 ... Batch8.HTO7 Batch8.HTO9
colData names(1): group_id
reducedDimNames(0):
altExpNames(0):

pb.logcounts = readRDS(file = "~/PFC_v3/PB_mean_logcounts_renormalized_subtypes.RDS")

meta.data = readRDS("~/PFC_v3/merged_metadata_with_SZTR.rds")

common.ids = intersect(colnames(pb.logcounts), meta.data$Internal_ID)
pb.logcounts = pb.logcounts[, common.ids]
metadata(pb.logcounts)$n_cells = metadata(pb.logcounts)$n_cells[, common.ids]

colData(pb.logcounts) = cbind(colData(pb.logcounts), meta.data)

Filter samples

ncells = metadata(pb.logcounts)$n_cells

mask = (fast_column_sums(ncells) > 100)
pb.logcounts.filtered = pb.logcounts[, mask]
metadata(pb.logcounts.filtered)$n_cells = metadata(pb.logcounts.filtered)$n_cells[, mask]

meta.filtered = meta.data[match(colnames(pb.logcounts.filtered), meta.data$ID), ]
colData(pb.logcounts.filtered) = cbind(colData(pb.logcounts.filtered), DataFrame(meta.filtered))

# Muscat analysis: two cohorts
require(muscat)
require(edgeR)
require(limma)

form = ~ Phenotype + Batch + Gender + Age + PMI + Benzodiazepines + Anticonvulsants + AntipsychTyp + AntipsychAtyp + Antidepress

resDE = lapply( levels(pb.logcounts.filtered$Cohort), function(chrt){

    keep.ids = colnames(pb.logcounts.filtered)[pb.logcounts.filtered$Cohort == chrt]

    pb.logcounts.filtered_sub = pb.logcounts.filtered[,keep.ids]
    metadata(pb.logcounts.filtered_sub)$n_cells = metadata(pb.logcounts.filtered_sub)$n_cells[,keep.ids]

    design.mat <- model.matrix(form, data = droplevels(colData(pb.logcounts.filtered_sub)))

    colnames(design.mat)[1] = c("Intercept")

    contrast.mat <- makeContrasts(contrasts = "PhenotypeSZ", levels = design.mat)

    pbDS(pb.logcounts.filtered_sub, method = "limma-trend", min_cells = 5, design = design.mat, contrast =  contrast.mat, filter = "gene")
})

2..6..16..19..21..25..30..2..6..16..19..21..25..30..

names(resDE) = levels(colData(pb.logcounts.filtered)$Cohort)

readr::write_rds(resDE, "~/PFC_v3/filtered_resDE_subtypes.rds")

celltypes = intersect(names(resDE$McLean$table$PhenotypeSZ), names(resDE$MtSinai$table$PhenotypeSZ))

filtered.tables = lapply(celltypes, function(celltype) {
  tbl1 = resDE[[1]]$table$PhenotypeSZ[[celltype]]
  tbl2 = resDE[[2]]$table$PhenotypeSZ[[celltype]]
  
  genes = intersect(tbl1$gene[1 <= abs(tbl1$logFC/sd(tbl1$logFC))], tbl2$gene[1 <= abs(tbl2$logFC / sd(tbl2$logFC))])
  
  tbl1 = tbl1[match(genes, tbl1$gene), ]
  tbl2 = tbl2[match(genes, tbl2$gene), ]
  
  mask = sign(tbl1$logFC)*sign(tbl2$logFC) > 0
  tbl1 = tbl1[mask, ]
  tbl2 = tbl2[mask, ]
  
  tbls = list(McClean = tbl1, MtSinai = tbl2) 
  tbls = lapply( tbls, function(tab){
    tab$se = tab$logFC / tab$t
    tab
  })
  
  return(tbls)
})
names(filtered.tables) = celltypes

readr::write_rds(filtered.tables, "~/PFC_v3/individual_diff_results_filtered_subtypes.rds")

Export as excel tables

for(ds in 1:2) {

  library(openxlsx)
  Up.wb <- createWorkbook()
  for(i in 1:length(tbls)) {
    res = filtered.tables[[i]][[ds]]
    res = res[res$logFC > 0, ]
    res = cbind(data.frame(Gene = rownames(res)), res)
    res = res[order(res$t, decreasing = T), ]
  
    n = names(filtered.tables)[[i]] #str_replace(arch.names[arch.order[i]], "/", "-")
    
    addWorksheet(wb=Up.wb, sheetName = n)
    writeData(Up.wb, sheet = n, res) 
  
  }
  saveWorkbook(Up.wb, sprintf("~/PFC_v3/DE_genes_up_%s_subtype.xlsx", names(filtered.tables[[i]])[[ds]]), overwrite = TRUE)
  
  library(openxlsx)
  Down.wb <- createWorkbook()
  for(i in 1:length(tbls)) {
    res = filtered.tables[[i]][[ds]]
    res = res[res$logFC < 0, ]
    res = cbind(data.frame(Gene = rownames(res)), res)
    res = res[order(res$t, decreasing = F), ]
    
    n = names(filtered.tables)[[i]] #str_replace(arch.names[arch.order[i]], "/", "-")
    
    addWorksheet(wb=Down.wb, sheetName = n)
    writeData(Down.wb, sheet = n, res) 
  
  }
  saveWorkbook(Down.wb, sprintf("~/PFC_v3/DE_genes_down_%s_subtype.xlsx", names(filtered.tables[[i]])[[ds]]), overwrite = TRUE)
}

combined.analysis.tables = lapply(names(filtered.tables), function(celltype) {
  print(celltype)
  tbls = filtered.tables[[celltype]]
  
  gene.tbls = lapply(1:nrow(tbls[[1]]), function(i) {
    dfs = lapply(1:length(tbls), function(k) tbls[[k]][i, ])
    df = do.call("rbind", dfs)
  })
  names(gene.tbls) = tbls[[1]]$gene
    
  combined.analysis.tbl = do.call(rbind, lapply(names(gene.tbls), function(gene){
    x = suppressWarnings(metafor::rma(yi=logFC, sei=se, data = gene.tbls[[gene]], method="FE"))
    combined.tbl = data.frame( gene = gene, 
        logFC     = x$beta,
        se        = x$se,
        tstat = x$zval,
        P.Value   = x$pval)
    return(combined.tbl)
  }))
  rownames(combined.analysis.tbl) = names(gene.tbls)
  
  combined.analysis.tbl = combined.analysis.tbl[order(combined.analysis.tbl$P.Value), ]
  
  return(combined.analysis.tbl)
})
names(combined.analysis.tables) = names(filtered.tables)

DF = do.call(rbind, combined.analysis.tables)
DF$adj.P.Val = p.adjust(DF$P.Value, "fdr")
combined.analysis.tables = split(DF, unlist(lapply(names(combined.analysis.tables), function(celltype) rep(celltype, nrow(combined.analysis.tables[[celltype]])))))

readr::write_rds(combined.analysis.tables, "~/PFC_v3/meta_analysis_diff_results_subtypes.rds")

resDE = readr::read_rds("~/PFC_v3/filtered_resDE.rds")

filtered.tables = readr::read_rds("~/PFC_v3/individual_diff_results_filtered.rds")

combined.analysis.tables = readr::read_rds("~/PFC_v3/meta_analysis_diff_results.rds")

  library(openxlsx)
  Up.wb <- createWorkbook()
  for(i in 1:length(combined.analysis.tables)) {
    res = combined.analysis.tables[[i]]
    res = res[(res$logFC > 0.1) & (res$P.Value <= 0.05), ]
    res = res[order(res$t, decreasing = T), ]
  
    n = names(combined.analysis.tables)[[i]] #str_replace(arch.names[arch.order[i]], "/", "-")
    
    addWorksheet(wb=Up.wb, sheetName = n)
    writeData(Up.wb, sheet = n, res) 
  
  }
  saveWorkbook(Up.wb, sprintf("~/PFC_v3/DE_genes_up_%s_subtype.xlsx", "combined"), overwrite = TRUE)
  
  
  library(openxlsx)
  Down.wb <- createWorkbook()
  for(i in 1:length(combined.analysis.tables)) {
    res = combined.analysis.tables[[i]]
    res = res[(res$logFC < -0.1) & (res$P.Value <= 0.05), ]
    res = res[order(res$t, decreasing = F), ]
    
    n = names(combined.analysis.tables)[[i]] #str_replace(arch.names[arch.order[i]], "/", "-")
    
    addWorksheet(wb=Down.wb, sheetName = n)
    writeData(Down.wb, sheet = n, res) 
  
  }
  saveWorkbook(Down.wb, sprintf("~/PFC_v3/DE_genes_down_%s_subtype.xlsx", "combined"), overwrite = TRUE)

DE.sc = matrix(0, nrow(pb.logcounts), length(combined.analysis.tables))
tstats = matrix(0, nrow(pb.logcounts), length(combined.analysis.tables))
logFC = matrix(0, nrow(pb.logcounts), length(combined.analysis.tables))
logPvals = matrix(0, nrow(pb.logcounts), length(combined.analysis.tables))
rownames(DE.sc) = rownames(tstats) = rownames(logFC) = rownames(logPvals) = rownames(pb.logcounts)
colnames(DE.sc) = colnames(tstats) = colnames(logFC) = colnames(logPvals) = names(combined.analysis.tables)

limma_trend_mean.scores = matrix(0, nrow(pb.logcounts), length(combined.analysis.tables))
Up.genes = vector("list", length(combined.analysis.tables))
Down.genes = vector("list", length(combined.analysis.tables))
rownames(limma_trend_mean.scores) = rownames(pb.logcounts)
names(Up.genes) = names(Down.genes) = colnames(limma_trend_mean.scores) = names(combined.analysis.tables)
for(i in 1:length(combined.analysis.tables)) {
    print(i)
    
    tbl = combined.analysis.tables[[i]]

    tstats[tbl$gene, names(combined.analysis.tables)[[i]]] = tbl$tstat
    logFC[tbl$gene, names(combined.analysis.tables)[[i]]] = tbl$logFC
    logPvals[tbl$gene, names(combined.analysis.tables)[[i]]] = -log10(tbl$adj.P.Val)
    
    DE.sc[tbl$gene, names(combined.analysis.tables)[[i]]] = tbl$tstat
    
}

[1] 1
[1] 2
[1] 3
[1] 4
[1] 5
[1] 6

limma_trend_mean.scores[is.na(limma_trend_mean.scores)] = 0


Up.genes = lapply(combined.analysis.tables, function(combined.analysis.tbl) {
  combined.analysis.tbl$gene[(combined.analysis.tbl$logFC > 0.1) & (combined.analysis.tbl$adj.P.Val < 0.05)]
})
Down.genes = lapply(combined.analysis.tables, function(combined.analysis.tbl) {
  combined.analysis.tbl$gene[(combined.analysis.tbl$logFC < -0.1) & (combined.analysis.tbl$adj.P.Val < 0.05)]
})

DE.new = list(DE.sc = DE.sc, tstats = tstats, logFC = logFC, logPvals = logPvals, Up.genes = Up.genes, Down.genes = Down.genes)
saveRDS(DE.new, "~/PFC_v3/DE_genes_pseudobulk_final_subtypes.rds")

ALS.DEs = read.table("~/als/MAGMA_DEG_input.tsv", sep = "\t")

suppressWarnings(ids <- AnnotationDbi::mapIds(org.Hs.eg.db, keys = ALS.DEs$V2, keytype = "ENSEMBL", column = "SYMBOL", multiVals = "first"))

'select()' returned 1:many mapping between keys and columns

ids[is.na(ids)] = ""
ids = as.character(ids)

ALS.genes = split(ids, ALS.DEs$V1)
X = do.call(cbind, lapply(ALS.genes, function(gs) as.numeric(rownames(variant_gene_scores) %in% gs)))
rownames(X) = rownames(variant_gene_scores)


# gg = apply(variant_gene_scores, 2, function(x) rownames(X)[scale(x) > 1])


EN = assess.geneset.enrichment.from.scores(variant_gene_scores, as(X, "sparseMatrix"))
logPvals = EN$logPvals
colnames(logPvals) = colnames(variant_gene_scores)


  require(ComplexHeatmap)
  pdf("~/als/MAGMA.pdf", height = 48)
Heatmap(logPvals)
  dev.off()

null device 
          1 

  require(ComplexHeatmap)
  pdf("~/als/MAGMA.pdf", height = 48)
  ComplexHeatmap::Heatmap(Res[, -5], rect_gp = gpar(col = "black"))

Warning in for (i in c(setdiff(seq_len(nc), c(1, nc)), c(1, nc))) { :
  closing unused connection 8 (~/als/magma/hmagmaAdultBrain__sz3_DE.gsa.out)
Warning in for (i in c(setdiff(seq_len(nc), c(1, nc)), c(1, nc))) { :
  closing unused connection 7 (~/als/magma/hmagmaAdultBrain__ms_DE.gsa.out)
Warning in for (i in c(setdiff(seq_len(nc), c(1, nc)), c(1, nc))) { :
  closing unused connection 6 (~/als/magma/hmagmaAdultBrain__alzKunkleNoapoe_DE.gsa.out)
Warning in for (i in c(setdiff(seq_len(nc), c(1, nc)), c(1, nc))) { :
  closing unused connection 4 (~/als/magma/hmagmaAdultBrain__alz2noapoe_DE.gsa.out)
Warning in for (i in c(setdiff(seq_len(nc), c(1, nc)), c(1, nc))) { :
  closing unused connection 3 (~/als/magma/hmagmaAdultBrain__als_DE.gsa.out)

  dev.off()

null device 
          1 

visualize.markers(In.ace, c("CALB1", "CALB2", "CUX2", "NOS1", "NPY"))

[1] "Markers Visualized: CALB1, CALB2, CUX2, NOS1, NPY"
[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

# Rosehip
# VIP
# REELIN
# 
# PV-Basker
# PV-Chan
# SST

epi_markers$Sst_Calb1

[1] "HPGD+"   "CALB1+"  "FBN2+"   "TAC3-"   "TH-"     "GXYLT2-" "ADGRG6-"

epi_annot = readr::read_rds("~/PFC_v3/epilepsy_paper_cell_annot.rds")
epi_markers = readr::read_rds("~/PFC_v3/epilepsy_paper_markers.rds")

In.cell.annot = annotate.cells.using.markers(In.ace, epi_markers[20:48])

Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5
Max it = 5

cc = grep("Sst", colnames(In.cell.annot$Enrichment))
sapply(cc, function(i) {
  plot(plot.ACTIONet.gradient(In.ace, x = In.cell.annot$Enrichment[, i], alpha_val = 0)+ggtitle(colnames(In.cell.annot$Enrichment)[[i]]))
})

            [,1]    [,2]    [,3]    [,4]    [,5]    [,6]    [,7]    [,8]    [,9]    [,10]  
data        List,0  List,0  List,0  List,0  List,0  List,0  List,0  List,0  List,0  List,0 
layers      List,1  List,1  List,1  List,1  List,1  List,1  List,1  List,1  List,1  List,1 
scales      ?       ?       ?       ?       ?       ?       ?       ?       ?       ?      
mapping     List,0  List,0  List,0  List,0  List,0  List,0  List,0  List,0  List,0  List,0 
theme       List,10 List,10 List,10 List,10 List,10 List,10 List,10 List,10 List,10 List,10
coordinates ?       ?       ?       ?       ?       ?       ?       ?       ?       ?      
facet       ?       ?       ?       ?       ?       ?       ?       ?       ?       ?      
plot_env    ?       ?       ?       ?       ?       ?       ?       ?       ?       ?      
labels      List,6  List,6  List,6  List,6  List,6  List,6  List,6  List,6  List,6  List,6 

annotate.profile.using.markers <- function (feature.scores, marker.genes, rand.sample.no = 1000) 
{
    if (is.matrix(marker.genes) | is.sparseMatrix(marker.genes)) {
        marker.genes = apply(marker.genes, 2, function(x) rownames(marker.genes)[x > 
            0])
    }
    specificity.panel = feature.scores
    GS.names = names(marker.genes)
    if (is.null(GS.names)) {
        GS.names = sapply(1:length(GS.names), function(i) sprintf("Celltype %s", 
            i))
    }
    markers.table = do.call(rbind, lapply(names(marker.genes), 
        function(celltype) {
            genes = marker.genes[[celltype]]
            if (length(genes) == 0) {
                err = sprintf("No markers left.\n")
                stop(err, call. = FALSE)
            }
            signed.count = sum(sapply(genes, function(gene) grepl("\\+$|-$", 
                gene)))
            is.signed = signed.count > 0
            if (!is.signed) {
                df = data.frame(Gene = genes, Direction = +1, 
                  Celltype = celltype, stringsAsFactors = FALSE)
            }
            else {
                pos.genes = (as.character(sapply(genes[grepl("+", 
                  genes, fixed = TRUE)], function(gene) stringr::str_replace(gene, 
                  stringr::fixed("+"), ""))))
                neg.genes = (as.character(sapply(genes[grepl("-", 
                  genes, fixed = TRUE)], function(gene) stringr::str_replace(gene, 
                  stringr::fixed("-"), ""))))
                df = data.frame(Gene = c(pos.genes, neg.genes), 
                  Direction = c(rep(+1, length(pos.genes)), rep(-1, 
                    length(neg.genes))), Celltype = celltype, 
                  stringsAsFactors = FALSE)
            }
        }))
    markers.table = markers.table[markers.table$Gene %in% rownames(specificity.panel), 
        ]
    if (dim(markers.table)[1] == 0) {
        err = sprintf("No markers left.\n")
        stop(err, call. = FALSE)
    }
    specificity.panel = specificity.panel[markers.table$Gene, ]
    IDX = split(1:dim(markers.table)[1], markers.table$Celltype)
    print("Computing significance scores")
    set.seed(0)
    Z = sapply(IDX, function(idx) {
        markers = (as.character(markers.table$Gene[idx]))
        directions = markers.table$Direction[idx]
        mask = markers %in% colnames(specificity.panel)
        A = as.matrix(specificity.panel[markers[mask], ])
        sgn = as.numeric(directions[mask])
        stat = A %*% sgn
        rand.stats = sapply(1:rand.sample.no, function(i) {
            rand.samples = sample.int(dim(specificity.panel)[2], 
                sum(mask))
            rand.A = as.matrix(specificity.panel[, rand.samples])
            rand.stat = rand.A %*% sgn
        })
        cell.zscores = as.numeric((stat - apply(rand.stats, 1, 
            mean))/apply(rand.stats, 1, sd))
        return(cell.zscores)
    })
    Z[is.na(Z)] = 0
    Labels = colnames(Z)[apply(Z, 1, which.max)]
    Labels.conf = apply(Z, 1, max)
    names(Labels) = rownames(specificity.panel)
    names(Labels.conf) = rownames(specificity.panel)
    rownames(Z) = rownames(specificity.panel)
    out = list(Label = Labels, Confidence = Labels.conf, Enrichment = Z)
    return(out)
}

names(curatedMarkers_human$Brain$PFC$Velmeshev2019$marker.genes)

 [1] "AST"              "Endothelial"      "IN-PV"            "IN-SST"           "IN-SV2C"          "IN-VIP"           "L2/3"            
 [8] "L4"               "L5/6"             "L5/6-CC"          "Microglia"        "Neu-NRGN"         "Oligodendrocytes" "OPC"             

Ast.DE = annotate.profile.using.markers(t(spec), mm)

[1] "Computing significance scores"

Heatmap(Ast.DE$Enrichment)